%% =========================== [Description] ================================
%% ---------------- implements the function -----------------
% Under multiple inputs
% Calculate the mean of the variance        
% Draw three variable output response curves in one graph             
% Draw response comparison curves for deterministic vs. random inputs
% ---------------- antecedent file -----------------
% quasi0StiffSuspensionCertain.slx
% quasi0StiffSuspension.slx
% n1nDpolyRoots.m
% GetStatistics.m/GetStat.m
% ------------------ enter -------------------
% Cylinder air pressure
% Piston distance from rear wall
% Cylinder node distance
% Effective area of cylinder
%% ====================== [Select excitation signal and output amount] ======================
clear
close all
disp('Pavement sine signal - input 1')
disp('Pavement step signal - input 2')
exc=input('Signal:');
disp('Output select acceleration - input 1')
disp('Output select dynamic deflection - input 2')
disp('Output select dynamic load - input 3')
outputAmount=input('outputAmount:');
swiMat=[-1,1];
excSwitch=swiMat(exc);
%% =========================== [parameter area] ==============================
d0Certain=32.5e-3; %% piston distance from back wall m
AcCertain=pi*0.1^2; % cylinder piston area m^2
L0Certain=0.467; % cylinder node distance m
kCertain=3; %cylinder air pressure atm
%--------------------------------------------------------------------------
mCertain=800;
m=800; %Mass on spring kg
m1=100; %Tire mass
Pa=101*10^3;
A_level=16; %-|
B_level=64; %-|-road
C_level=256; %-|
Level=B_level;
v=60; %Vehicle speed
C=2000; %damping
qA=0.1; % excitation frequency
f=1;
%--------------------------------------------------------------------------
D=1; %number of random variables
coeOfVarCell={[0.1,0,0,0],[0,0,0.1,0],[0,0,0,0.1],...
    [0.2,0,0,0,0],[0,0,0.2,0],[0,0,0,0.2]};
MCnumber=500;
vMat=[d0Certain,AcCertain,L0Certain,kCertain];
order=3; 
%% ========================= [Calculation area] ================================

%% ============================[C] ==================================

sim('quasi0StiffSuspensionCertain.slx')
Ctable={accCertain,devCertain,LoadCertain};
opDataCertain=Ctable{outputAmount};

%% ===========================[PCE] =================================
pointNumber=200;
PCmean=zeros(pointNumber,6);
PCvar=zeros(pointNumber,6);
Collocation=n1nDpolyRoots(order,D); 
groupNumber=length(Collocation(:,1));
for ii=1:1
    coeOfVar=coeOfVarCell{ii}; 
    for groNum=0:groupNumber  
        if groNum==0
            [d0,Ac,L0,k]=iniRanVar(d0Certain,AcCertain,L0Certain,kCertain,...
                         coeOfVar,Collocation,groupNumber,1);
            sim('quasi0StiffSuspension.slx')
            outputData=zeros(length(a.time),groupNumber);
        else
            [d0,Ac,L0,k]=iniRanVar(d0Certain,AcCertain,L0Certain,kCertain,...
                         coeOfVar,Collocation,groupNumber,groNum);
            sim('quasi0StiffSuspension.slx')
            switch outputAmount %Select the output amount
                case 1 
                    output=a.data;
                     yLabelText='\it{Acceleration}$\rm{(m/s^2)}$';
                    proUnit='$\rm{(s^2/m)}$';
                case 2
                    output=fd.data;
                    yLabelText='\it{Deflection}$\rm{(m)}$';
                   proUnit='$\rm{(1/m)}$';
                case 3
                    output=Fd.data;
                    yLabelText='\it{Relative Dynamic Load}';
                    proUnit='';
                otherwise
                    error('Select 1\2\3');
            end
            outputData(:,groNum)=output;
        end
       %-----------------------------------------------------------------------
        proPers=floor(groNum/groupNumber*1); %
        clc %---- progress bar area
        disp(['PCE[','#'*ones(1,proPers),' '*ones(1,100-proPers) ...          %
            ,']',num2str(proPers),'%']) %
       % -----------------------------------------------------------------------
    end
    timeNumber=10;
    sampleTime=length(a.time)-1;
    stepL=timeNumber/pointNumber;
    time=0:stepL:timeNumber-stepL; 
    for i=1:pointNumber    
        [PCmean(i,ii),PCvar(i,ii)]=GetStat ...
        (Collocation,outputData((sampleTime/pointNumber*(i-1))+1,:)',order); 
        %----------------------------------------------------------------------
        proPers=floor(i/pointNumber*99); %
        clc %---- progress bar area
        disp(['PCE[','#'*ones(1,1),'#'*ones(1,proPers),' '*...                %
            ones(1,99-proPers) ,']',num2str(1+proPers),'%']) %
        % ----------------------------------------------------------------------
    end
end
%% ==================== [Chaotic mean and model mean plot lines] ======================
close all;
edgeCell={'h-','s-','d-','o-','^-','v-'};
colormap=[ 0 0.4470 0.7410;
    0.8500 0.3250 0.0980;
    0.9290 0.6940 0.1250;
    0.4940 0.1840 0.5560;
    0.4660 0.6740 0.1880;
    0.3010 0.7450 0.9330;
    0.6350 0.0780 0.1840];
figure(2)
plot(opDataCertain.time(1:25:10001),opDataCertain.data(1:25:10001),':','LineWidth',1.5,'color','red')
hold on
for jj=1:6
    plot(time,PCmean(:,jj),edgeCell{jj},'color',colormap(jj+1,:))
    hold on
end
if exc==2
    xlim([0,5])
end
ylabel(yLabelText,'Interpreter','latex','FontSize',20)
xlabel('Time(s)','Interpreter','latex','FontSize',20)
set(gca, 'Fontname', 'Times New Roman','FontSize',20)
legend('$\zeta=0$','$d_0(\zeta=0.1)$','$L_0(\zeta=0.1)$','$P_{c0}(\zeta=0.1)$','$d_0(\zeta=0.2)$','$L_0(\zeta=0.2)$','$P_{c0}(\zeta=0.2)$','$P_{c0}(\ zeta=0.2)$','Interpreter','latex')

%% ========================== [mean of variance] =============================
PCvarmean=mean(PCvar,1);
disp(PCvarmean)
PCmeanmax=max(PCmean);
disp(PCmeanmax)
Certainmax=max(opDataCertain.data);
disp(Certainmax)

